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ABSTRACT 

Mutual sticking of dust aggregates is the first step toward planetesimal formation in protoplanetary disks. In 
spite that the electric charging of dust particles is well recognized in some contexts, it has been largely ignored 
in the current modeling of dust coagulation. In this study, we present a general analysis of the dust charge state 
in protoplanetary disks, and then demonstrate how the electric charging could dramatically change the currently 
accepted scenario of dust coagulation. First, we describe a new semianalytical method to calculate the dust 
charge state and gas ionization state self-consistently. This method is far more efficient than previous numerical 
methods, and provides a general and clear description of the charge state of gas-dust mixture. Second, we 
apply this analysis to compute the collisional cross section of growing aggregates taking their charging into 
account. As an illustrative example, we focus on early evolutionary stages where the dust has been thought 
to grow into fractal (D ~ 2) aggregates with a quasi-monodisperse (i.e., narrow) size distribution. We find 
that, for a wide range of model parameters, the fractal growth is strongly inhibited by the electric repulsion 
between colliding aggregates and eventually "freezes out" on its way to the subsequent growth stage involving 
collisional compression. Strong disk turbulence would help the aggregates to overcome this growth barrier, 
but then it would cause catastrophic collisional fragmentation in later growth stages. These facts suggest that 
the combination of electric repulsion and collisional fragmentation would impose a serious limitation on dust 
growth in protoplanetary disks. We propose a possible scenario of dust evolution after the freeze-out. Finally, 
we point out that the fractal growth of dust aggregates tends to maintain a low ionization degree and, as a result, 
a large magnetorotationally stable region in the disk. 

Subject headings: dust, extinction — methods: analytical — planetary systems: formation — planetary sys- 
tems: protoplanetary disks — plasmas 
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1. INTRODUCTION 

The initial step toward planetesimal formation in protoplan- 
etary disks is the collisional growth of submicron dust grains 
into macroscopic aggregates. A standard scenario is that dust 
aggregates grow by mutual sticking, gradually settle to the 
midplane of the disk, and finally form a dense dust layer. It is 
still an open issue whether subsequent growth is established 
by the gravitational instability of the layer or the direct growth 
of the aggregates. To address this issue, further understanding 
on earlier evolutionary stages is needed. 

It has been recognized that the internal structure of ag- 
gregates is a key factor for their growth and settling. Early 
studies on dust coagulation model ed the aggregate s as a 
compact, n onporous object (e.g., Weidenschi llingI 119801: 
iNakagawae t al. 1981). Both numerical simulations and lab- 
oratory experiments have revealed, however, that aggregates 
are not at all compac t, bu t has an open, fluffy structure (for a 
review, see Meakin 1 9911 Ilium 2004; Dominik et al1l2007h . 
This is particularly true for aggregates formed at an early 
growth stage where the collisional velocity is too low for col- 
liding aggregates to compress each othe r. It has been ob- 
served in numerical (iKempf et al. ill999l) as weU as experi- 
mental (IWurm & Blu m 1998; Blu m et al.lll998ll2000l) studies 
that the outcome is an ensemble of fractal aggregates with the 
fractal dimension D <2 and with a quasi-monodisperse (i.e., 
narrow) mass distribution. This fractal growth typically lasts 
until the aggregates become centimeter-sized (ISuvama et al.l 
|2008|) . A remarkable dynamical property of these fluffy ag- 
gregates is that they keep a strong couphng to ambient gas and 
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thus a low drift velocity relative to the gas throughout the evo- 
lution. This could be crucial to the formation of very thin dust 
layer where planetesimals may be formed by gravitational in- 
stability. 

Dust grains and aggregates are not only the building block 
of planetesimals but also powerful absorbers of charged par- 
ticles in the gas disks. It is now widely accepted that tur- 
bulence in th e disks is attributed to m agnetorotational insta- 
bility (MRI; iBalbus & Hawlevl[T99ll) . For this mechanism 
to work, at least a part of the disk needs to be sufficiently 
ionized for the gas to couple to magnetic fields. Many au- 
thors have examined whether protoplanetary disks can be 
ionized enough to sus t ain MHD turbulen ce (Gammie 199^ 
Sass2old et al.' '19971; 'Sano et al.' '20001; llgea & Glassgoldl 
1999;'llgner & Nelson 2006a b e; Wardii|2007|). One of the 
important findings is that the turbulent region is strongly con- 
trolled by the concentration of dust materials since they ef- 
ficiently remove away ionized particles from ambien t gas 
(ISano et al.ll2000l; lllgner & Nelsonll2006al IWardldllOOTh . 

Although the importance of dust charging is well recog- 
nized in the above context, its effect on dust coagulation 
in protoplanetary disks has been hardly examined. Charg- 
ing of aggregates causes electrostatic interaction between 
them, which may significantly increase or decrease the co- 
agulation rates. Recently, a series of studies have suggested 
that charge-induced dipole interaction might trigger runaway 
growt h of dust aggregates (Ivlev et al. 2002; Konopkaet al] 
l2005l) . However, these studies considered a situation where 
the ambient gas is not ionized and the net charge of dust ag- 
gregates vanishes identically. In protoplanetary disks, on the 
contrary, the net charge of dust aggregates does not vanish 
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due to the presence of weakly ionized ambient gas, and there- 
fore both dust charging and gas ionization must be taken into 
account. 

This study explores how the electrostatic charging of dust 
aggregates could be crucial to their coagulation in protoplan- 
etary disks. For this purpose, we have to know in advance 
how the charge state of aggregates evolves with their growth. 
This is a complicated problem, since we also have to solve 
the ionizati on state of ambient gases self-consistently. Previ- 
ous st udies dSano et al.ll2000l:lllgner & Nelson 20Q6a; Wardld 
l2007h have handled this problem with direct numerical cal- 
culations in which dust particles with different charges and 
sizes are treated as different charged species as well as many 
species of ions. However, this approach becomes inefficient 
when one tries to solve this problem and dust growth simul- 
taneously, since the dispersion of charge and size increases as 
the dust growth. The central strategy taken in this study is to 
solve this problem as analytically as possible. This approach 
does not only reduce the computational expense but also pro- 
vides general insight into the charge state of gas-dust mixture. 
As a result, we show that all the conditions for ionization- 
recombination equilibrium are reduced to a single algebraic 
equation. Just by solving this equation numerically, we can 
obtain both of the dust charge state and the gas ionization state 
analytically. We also confirm that the semianalytical calcula- 
tions agree very well with direct numerical calculations using 
the original equations. This semianalytical method will be a 
powerful tool for the simulations of charged dust coagulation 
and MRI turbulence. 

As an illustrative example, we calculate the collisional 
cross section of dust aggregates growing in a protoplane- 
tary disk taking into account their electric charging. We fo- 
cus on early stages of dust evolution where the aggregates 
has been tho ught to experience fractal, quas i-monodisperse 
growth (e.g.. iBlumlllOM: iDominik et al.ll2007l) . For a wide 
range of model parameters, we find that the effective cross 
section is quickly suppressed as the fractal growth proceeds 
and finally vanishes at a surprisingly early stage. This means 
that the fractal growth "freezes out" on its way to the sub- 
sequent growth stage where collisional compression of ag- 
gregates occurs. This is because the electrostatic repulsion 
between aggregates becomes strong enough to prevent their 
mutual collision. Strong turbulence in the disk will help the 
aggregates to overcome this electric barrier, but then it will 
cause catastrophic disruption of collided aggregates at later 
stages. Therefore, if the freeze-out of the fractal growth truly 
means the end of dust evolution, the combination of the elec- 
tric charging and the collisional disruption imposes a very 
strict limitation on dust coagulation and subsequent planetes- 
imal formation in protoplanetary disks. Our findings strongly 
suggest that the dust charing effect should be seriously taking 
into account in the modeling of dust evolution. 

This paper is organized as follows. In §2, we present a set 
of equations that describes the reactions of charged particles 
(ions, electrons, and dust aggregates), and derive the equa- 
tion that determines the equilibrium state. In §3, we calculate 
the electrostatic repulsion energy between two colliding ag- 
gregates to show that the quasi-monodisperse fractal growth 
is strongly inhibited for a wide range of disk parameters. In 
§4, we discuss the validity of some important assumptions 
and point out a possible scenario of dust evolution after the 
"freeze-out" of the fractal growth. A summary is presented in 
§5. 
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Fig. 1 . — Schematic diagram of tlie ionization-recombination reactions in a 
gas-dust mixture. Ions and free electrons are created (solid arrows) by some 
ionizing sources (e.g., cosmic rays) and are removed tlirougli the gas-phase 
recombination (dotted arrows) or the adsorption to dust (double line arrows). 
Some species of ions may react with neutral gas particles to create different 
species of ions (dashed arrows). The equilibrium charge distribution ni(Z) 
of dust aggregates ai'e determined by the balance between all these reactions. 



2. EQUILIBRIUM CHARGE DISTRIBUTION 

2.1. Kinetic equations for ionization-recombination 
reactions 

We model the ionization-recombination reactions in a gas- 
dust mixture as follows (see also figUli- Some ionizing 
sources (e.g. cosmic rays) create ions X|*' and free electrons 
from neutral particles X^*' at a rate C^*^* (here k(= 1,2, • • •) la- 
bels each species of ions and associated neutrals). We assume 
that the ions and electrons are quickly thermahzed and have 
thermal velocities uf^ and Wg, respectively. We neglect the 
possibility that the free electro ns might be much mor e ener- 
getic in a MRI-active region (llnutsuka & Sanoll2005h . The 
ions may react with neutrals Xg* '' to produce another species 
of ions X^ (Z ^k), or may recombine with free electrons in the 
gas phase. We denote the rate coefficient for the ion-neutral 
reaction and the gas-phase recombination by /3 ^'^''^ and 
respectively. Also, the ions and free electrons may collide 
with dust aggregates to adsorb onto their surfaces. We write 
the collisional cross section for an aggregate and an ion (elec- 
tron) as CTdi(de) ■ These cross sections generally include the ef- 
fect of electrostatic interaction as well as the sticking proba- 
bility (see §2.2). Each dust aggregate may have different in- 
ternal structure and charge Ze from the other. We represent a 
set of parameters describing the structure (e.g., mass, radius) 
as / = {/i,/2, • • •}. In this section, we assume that the above 
reactions proceeds faster than the mutual collision of dust ag- 
gregates and treat / as constant parameters. The validity of 
this assumption is discussed in §4.1. 

The above charge reactions are described by a set of kinetic 
equations. Let us denote the number densities of ions X|*', 
free electrons, and dust aggregates as nf\ n^, and ni(I,Z), 

respectively. The rate equations for nf\ n^, and nd(I,Z) are 
given by 

«f = CW«W-«W«« / dlJ2a,,(I,Z)na(I,Z) 

z 
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and 

nd(/,Z) = ^MW«fVdi(/,Z- l)nd(/,Z- l)-adi(/,Z)«d(/,Z)] 

+Mene [crde(/, Z+ 1 )nd(/, Z + 1 ) - (Tde(/, Z)nd(/, Z)] , (3) 

respectively. In equations ([TJ and (l2]i, n^*^' and Mg*^ '' denote 

the number densities of neutrals and xf ''. We assume 
that neutral particles are much more abundant than charged 
particles and regard nS^^ and n!f as constant parameters. 
In addition, we impose the charge neutrality condition 



^nf^-ne+ / <i/^Z«d(/,Z) = 0. 



(4) 



Equations ([T]i-(|4|i form the closed set of basic equations for 
ionization-recombination reactions in gas-dust mixture. 

2.2. The equilibrium solution 

The equilibrium solution is obtained by imposing the con- 
ditions 

(5) 



nf' = rie = rid(I,Z) = 

for all k, I, and Z. 

Usually, the equilibrium solutions are calculated with a re- 
action schem e involvin g many i on species and ion-neutral re- 
actions (e.g., lUmebayash i & Nakanol ll980l: lllgner & NelsonI 
l2006a) . Because of the complexity of the ion-neutral reac- 
tions, it is generally impossible to solve this problem analyti- 
cally without any approximation or simplification. 

This problem, however, becomes analytically tractable if 
we dot not require to distinguish ion species. This is just 
achieved by introducing the total ion density 



(6) 



Taking the sum of equations ([T]l over all k, we obtain the rate 
equation for rii, 



(. It is evident that equation (I?) is much simpler than the orig- 
inal equation ([U. This is mainly attributed to the cancellation 
of the last term in the original equation, i.e., the term describ- 
ing the ion-neutral reactions. Of course, this also means that 
we have lost the chance to know the composition of ions in 
detail. This fact does not bother us since our primary interest 
is the charge state of dust aggregates, not the composition of 
ions. 

Now we try to solve the equations (|2|i-(|4|i and (|7]i under the 
equilibrium condition (|5]l as analytically as possible. First, 
equations dTji and (|2]i are written as 



CWg - Mi (cTdi ) nd«i - /^MiHe = , 



Cng-Me((Tde)ndne-/3nine = 0, 



(11) 

(12) 



respectively. Here we have defined the averages of an arbi- 
trary function F = F{I,Z) over Z and over / as 



{F)il)-: 



F(Z) = 



1 



nd(/) Y 



F(I,Z)na{I,Z) 



F(I,Z)nA(I,Z)dI 



(13) 



(14) 



with tiiil) = 'Y^^Hdil ,Z) and nd(Z) = j nd{I,Z)dI, respec- 
tively. We have also defined the total number density «d of 
dust aggregates by = "Y^-^j ^<i(Ji2)dL Equations (flTT i and 
(fT2l i can be easily solved in terms of ni and «e as 



Me(0-de)«d 
2/3 

Mi(o-di)nd 



4/3Cng 



MiMelCTdi) (Crde)nd 



4/3C«g 



■1 



(15) 



(16) 



2/3 \ V 
Next, equation ^ is reduced to 

MiniO-di(/,Z- l)nd(/,Z- l)-MeneCTde(/,Z)nd(/,Z) 

= MiniCrdi(/,Z)nd(/,Z)-MeneO-de(^,Z+ l)nd(AZ+ 1). (17) 

This equation means that the "flux" of the distribution tiiil^Z) 
from the state Z - 1 to Z must be balanced to that from Z to 
Z + 1. If the flux did not vanish, there would exist a steady 
"flow" of the charge state distribution nd(I,Z) streaming from 
one direction to the other in Z-space. However, since nd{I,Z) 
must be vanish at Z ^ ±oo, such a steady flow must not ex- 
ist. Therefore, both of the left and right-hand sides of equa- 
tion {Tt\ must be zero. Hence we have 



«i = Cng-Mini / dlYadi{I,Z)nd{I,Z)-(3nin^, (7) UiniadiiI,Z)na{I,Z) = u^n^cra^iI,Z+l)nA{I,Z+l) 

J -y 



(18) 



where «g = J2k '■^'■^l number density of the gas, and 



1 (k) I 

Ui = — y n\ 'u. 



(8) 
(9) 
(10) 



are the average ion velocity, gas-phase recombination rate co- 
efficient, and ionization rate, respectively. Similarly, equa- 
tions ©-dill can be also written down in terms of n^, «[, /3, and 



This is just the condition of detailed balance between charge 
states Z and Z+ 1. 

For the sake of later convenience, we here rewrite the 
charge neutrality condition (|4| using the definitions (|6]l, ( fT3] ), 
and ( fT4b as 

«i-ne+"(zynd = 0. (19) 

As we will see later, this is the final equation that determines 
the equilibrium solution. 

The next step is to solve the detailed balance equation dTsl ). 
To carry out the calculation, we need to specify the forms of 
the effective collision cross sections, (Jdi and a^^. For simplic- 
ity, we model a fractal aggregate as a spherical, porous body 
with radius a and projected cross section a. Also, we ne- 
glect the electric polarization of aggregates (iDraine & SutinI 
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[19871) since a fractal aggrega te is likely to have low dielec- 



tricity. iDraine & SutinI (Il987h found that the strength of the 



electric polarization relative to the charge-charge (Coulomb) 
interaction is determined by a factor (e- l)/(e + 2), where e 
is the dielectric constant of the material. Given that a porous 
aggregate is well approximated by a dilute medium, the fac- 
tor (e- l)/(e+2) is inversely proportional to i ts mean interna l 
density (Clausius-Mossotti relation; see, e.g.. lJacksonll 19981) . 
For example, the mean internal density of a fractal aggre- 
gate with D w 2 is about inversely proportional to its radius 
a. Therefore, we expect that the polarization effect of such a 
low-density aggregate is negligible. Neglecting the polariza- 
tion effect, th e effective cross sections are simply written as 
(ISpitzeiill941h 



crdi(/,Z) = cr«i(Z) X 



crde(/,Z) = crSe(Z)x 




Z<0, 



Z>0, 



Z<0, 



z > 0, 



(20) 



(21) 



where >Si(e)(Z) is the probability that a colliding ion (electron) 
sticks to one of constituent monomers, and A = /k^T . In 
this study, we assume that s-^Z) and Si.(Z) are independent of 
the net charge Z carried by an aggregate, i.e., ii(e)(Z) = s^^). 

As shown in Appendix, the solution 11^(1, Z) to equation 
(flSl ) with equations ( l20l l and (l2Tl i is well approximated by a 
Gaussian distribution 



«d(/,Z): 



nd(/) 



v/2^(AZ2)„ 



exp 



(Z- (Z) J2 



2(AZ2)„ 



where 



nd(/) 



^Z«d(/,Z); 



and 



(AZ^ 



nd(/) 



^(z-(z)„)V/,z) = 



Ta 
'T 

l+T a 



2 + r A 



(22) 



(23) 



(24) 



are the mean and dispersion of the charge distribution for 
fixed a, respectively. This solution is valid when the radius 
a is much larger than A, as is for aggregates much larger than 
constituent monomers (see Appendix). The nondimensional 
parameter r = -(Z) A/a = {-{Z)e^ /a)/kBT measures the elec- 
trostatic attraction (repulsion) energy between a charged ag- 
gregate and an incident ion (electron) relative to the thermal 
kinetic energy. 

Since the charge distribution is parametrized by F only, Wi 
and He can be written as a function of a single parameter F. 
Using equations (I22l l- (l24l l. {an) and (crde) are evaluated as 

:CT,?,(l+r), (25) 



((Tde ; = CTie exp 



-F + 



A(l + F) 



2fl(2 + F) 



:CT5eexp(-F), (26) 



respectively. Here we have used in equation ( l26l l that a/ A ^ 
1. Substituting equations ( |25] | and (|26] | into equations ( fTSl l 
and ( fTSI l. we have 



«i = 



Cng Vl+2g(F)-l 
SiUiarii {\ + V)g{V) 



(27) 



where 



^(r) = 



C«g Vl+2g(F)-l 

ieMeCTHd eXp(-F)g(F) 

2l3(^ng expF 



(28) 



(29) 



SiUiSeUgian^y l+T 

Finally, the neutrality condition ( fT9l ) with equations ( |23]) . 
( |27] |. and ( |28] | leads to the equation for F, 



1 



l+F 



SiUi 1 

expF+- 



■0, 



e vi+2g(F)-i. 

where we have defined a nondimensional parameter 



6; 



(30) 



(31) 



We have consequently arrived at the conclusion that all the 
conditions for ionization equilibrium, equations (fTTT i. (fT2l) . 
and ( fTSl) . are reduced to a single equation ( l30l l for a single 
parameter F. 

To summarize the above analysis, we have considered the 
charge state of a dust-gas mixture in the presence of ionization 
sources. We have found that the self-consistent equilibrium 
solutions are written as analytical functions of a single param- 
eter F (eqs. ||22l-||24l, 1271, and ||28]), and have obtained the 
equation for this master parameter (eq. 1301 ). This equation 
can be easily solved numerically, and thus gives semianalyti- 
cal solutions to tii, n^, and nd{I,Z). 

The resultant equations can be further simplified when the 
gas-phase recombination rate coefficient f3 is so small that the 
factor g defined in equation (|29] | is much less than unity. In 
this limit, equations ( |27] |. (|28] ). and (l30l ) are simply rewritten 
as 

Cng 1 



SiMiCrnd l+F' 



:^expF, 



and 



l + F 



SiMi F 

expF+- 



^0, 



(32) 



(33) 



(34) 



respectively. As seen in §3.1.3, this approximation is valid 
in typical protoplanetary disks unless the dust is significantly 
depleted (e.g., by vertical sedimentation). 

2.3. Limiting cases 

Equation (|30] l. or ( [34] i. provides clear insight into the charge 
state of a gas-dust mixture. The first and second terms in the 
bracket in this equation originate from rig and -{Z)ni in the 
charge neutrality condition (01, respectively. This means that 
the parameter determines which of free electrons and dust 
aggregates are the dominant carriers of negative charge. In the 
following, we categorize the charge state from limiting cases 
of equation ( l34b . 

2.3.1. O oo: ion-electron plasma limit 

In the limit 8 oo, the contribution from charged dust 
becomes negligibly small. Hence, equation ( l34b is well ap- 
proximated by 

' ^'"^ - (35) 



l+F 



il^expF. 

SpUp 



This is a well-known charge-equilibrium condition fo r a dust 
particle immersed in an ordinary ion-electron plasma dSpitzeii 
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119411; IShukla & MamuDll2002l) . We denote the solution T to 
equation ([35]) by Fmax, since this is the maximum value of F 
obtained from equation (|34] |. Typically, Fmax takes a value of 
order unity. 

2.3.2. —>■ 0: ion-dust plasma limit 

In the limit 9 ^ 0, on the other hand, the contribution from 
free electrons becomes negligible. This means that the dom- 
inant carriers of negative charges are dust particles, not free 
electrons. We shall refer to this limit as the ion-dust plasma 
limit. In this hmit, equation (|34] | can be approximated by 

The solution is easily obtained as 



F; 



Vl +402-1 



e. 



(37) 




Since F is now negligibly small, the effective cross sections 
(cdi(e)) are approximately equal to asi^e). Equations (|32] | and 
( |33] ) lead to the ratio of n^ to ni, 



(38) 



where we have used that Mi(e) = y/^k^ T / nmn^) . We note t hat 
this value of We/ni is larger than that of lUmebayashil (1 1 983h by 
a factor of 1 / ^fs^. This is becaus e we have assumed t he value 
of ie as independent of Z while lUmebavashil ( Il983h consid- 
ered ie(Z > 0) = 1. As far as the author knows, there is no ex- 
perimental data that validates either of the assumptions. How- 
ever, this difference is practically unimportant unless the ratio 

Sg{Z < 0) /s^iZ > 0) is much less than unity. 

The value of O at the transition from one plasma regime to 
the other can be roughly estimated by equating the asymptotic 



solutions for both limits, i.e., 8 
estimation in the next section. 



We will use this 



3. APPLICATION: ELECTRIC BARRIER AGAINST 
DUST GROWTH 

Electrostatic interaction between charged aggregates affect 
their collisional cross section. Let us consider two dust aggre- 
gates with mass mj, radius aj and charge Zje, where j{= 1 , 2) 
labels the aggregates. The kinetic energy for relative motion 
of two aggregates 1 and 2 is written as 



1 2 



(39) 



where m = m\m2/(m\ +OT2) and Am are the reduced mass and 
the relative speed for the aggregate pair. The electrostatic en- 
ergy between the aggregates just before contact is 



^el ■ 



(40) 



CTdd = 



a\ +fl2 

Neglecting the polarization effect as done in §2, the effective 
collision cross section (Tdd for the aggregates is expressed as 
(ISpitzeiill941l) 

7r(ai+a2)^ ^1--^^ , -Ekin > -Eel, ^4^^ 

0, E'kin < -Eel, 

Therefore, the condition for the aggregates to colhde with 
each other is 

£kin>-Eel. (42) 

In this section, we examine whether this condition is satisfied 
in an early stage of dust evolution in a protoplanetary disk. 



3.1. Model setup 

3.1.1. Protoplanetary disk model 

We assume that the gas surface density Eg and the temper- 
ature T of the disk obey power laws 



Eg(r)=1.7xlOVE 



lAU 



-3/2 



g/cm^, 



and 



T{r) = 280 



(lAu) 



-1/2 



K, 



(43) 



(44) 



where r is the distance from the central star and /s is a nondi- 
mensional scaling parameter The model with /s = 1 is known 
as the minimum-mass solar nebula (MMSN) model (Hayash^ 
Il98ll) . We adopt /s = 1 unless otherwise noted. The tem- 
perature profile (|44] | is valid only for optically thin regions. 
We nevertheless employ this profile throughout the disk since 
the main result is insensitive to the detail of the temperature 
profile (see eq. ll65l below). 

The hydrostatic equilibrium in the vertical direction gives 
the gas density distribution 



/2TrH(r) 



exp 



2H{ry 



(45) 



where z is the height from the disk midplane and H(r) = 
Cs(r)/il.Kir) is the gas scale height. The isothermal sound ve- 
locity Cj and the Kepler rotational frequency ilKir) are given 
by Cs(r) = ^k^Tir)/ pmii and i^Kir) = y/GM^/P, where p, is 
the mean molecular weight, niu is the hydrogen mass, G is the 
gravitational constant, and is the mass of the central star. 
We adopt fj, = 2.34 and M» = IMq in the following calcula- 
tion. The total number density of gas particles is given by 

ng = Pg//ZOTH. 

We assume that dust material is well mixed in the disk and 
that the dust density is simply related to the gas density pg 
by 

Pd(r,z) = fdgPg(r,z), (46) 

where /dg is the dust-to-gas ratio in the disk. The solar sys- 
te m abundance of cond ensates including water ice estimated 
by iPollack et"an (1 19941) leads to the dust-to-gas ratio /dg = 
0.014 as well as th e monomer bulk density po = 1.4g/cm^ 
dTanaka et al.ll2005l) . We do not consider the sublimation of 
water ice in inner disk regions for simplicity. 

3.1.2. Ionization rate 

In this study, we co nsi der Galactic cosmic 
rays (Umebavash i & Nakanol Il981h . stellar X- 
rays (Igea & Glas sgold 



19991) . and radionuclides 
dUmebayashi & Nakanol f2009r as the ionizing sources. 



(47) 



Thus, we decompose the ionization rate C as 

C ~ CCR + CxR + CrA 

where (qr, Cxr, and Cra denote the rate of ionization by cos- 
mic rays. X-rays, and radionuclides, respectively. We do not 
consider thermal i onization. This is ne gligible for T <C lO^^K, 
or for r '3> 0.1 AU (lUmebavashill 1 9 8 3l) . Charged particles are 
created primarily by ionization of H2 and He. The ioniza- 
tion rate for He is related to that for H2 by C'"''' = 0.84C<"=* 
(Umebavashi & Nakano 1990, 200^, so it is sufficient to 
know C*^^' only. The total ionization rate ( is given hy ( = 
C'^^-'jch, + C^^'^^'^He, where X112 = nw^/n^ and xhe = nue/ng are 
the fractional abundances of H2 and He. We calculate jch, and 
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xue fro m the solar system abundance by lAnders & Grevess"^ 
([T989h . 



The cosmic -ray ionization rate Crp-' for H2 is given by a 



fitting formula (Umebavashi & Nakano 2009) 



AHl) 
SCR.O 



exp 



XCR 



1 + 



XCR 



+ exp I 



Xg(r,z) 

XCR 



1 + 



x;('-,2)\3/4 



/Xgl^\ 
^ XCR ^ 



-4/3' 



-4/3 



,(48) 



where Ccro ~ 1-0 x 10 '^/s is the cosmic-ray ionization rate 
for H2 in the interstellar space, xcr ~ 96g/cm^ is the attenua- 
tion length of the ionization rate, and xti'',z) = /, Pg{r,z')dz' 

and x^('",-z) = ^g('") ~ Xg('":-z) ^re the vertical gas column 
densities measured from the upper and lower infinities, re- 
spectively. For the radionuclide ionization rate, we assume 

c, 



(H2) 

RA 



7 X 10 /s, which corresponds to the ionization 
rate by ^^Al with an ab undance ratio -*A1/^^A1 = 5 x 10"^ 
(lUmebayashi & Nakanoll2 009). 

The stellar X-ray ionization rate has been calculated by 
Hgea & Glassgold (1999) using the Monte Carlo radiative 
transfer code including Compton scatter ing. A useful fitting 
formula is given by Turner & Sand (l2008i) . 

CxR('",2)~CxR.of 7777) (7, 77^^ r) 

VIAU/ \2 X lO-^^erg/s/ 



exp 



lAUy V2xl030erg/s 

Xg('-,z)\ / Xs('",2) 



XXR 



-I- exp 



,(49) 



XXR 

where Lxr is the X-ray luminosity, and Cxr,o = 2.6 x 10"'^/s 
and xxr = 8.0g/cm^ are the fitting parameters. This fitting 
fo rmula approximately rep roduces the A;b7xr = 5keV result 
of llgea & Glassgoldl (119991) within the column density range 
Xg ^ Ig/cm^, where scattered hard (> 5keV) X-rays are re- 
sponsible for the ionization. We use equation (|49] | in the 
following calculation since the typical value of Xg is within 
the above range. We take Lxr = 2 x lO-'^erg/s in accordance 
with the median characteristic X-ray luminosity observed 
by Chandra for young solar-mass stars in the Orion Neb- 
ula Cluster (IWolk et al.l 12005). Although the characteristic 
X-ray temperature A;b7xr ~ 2.4keV observed bv IWolk et al.l 
(l2005h is lower than the assumed value of k^TxR « 5keV, 
the choice of the temperatur e does not significantly affect 
the resulting ionization rate (ilgea & Glassgoldl Il999l) . We 
do not consider temp orary increase in the X-ray luminosity 
due to stellar flaring (IWolk et al.l l2005). since our analytical 
metho d assumes station a ry ioniz ation processes. As pointed 
out by lllgner & NelsonI (l2006d) . the flaring could quantita- 
tively change the ionization state of the disk. We will examine 
the effect of the time-dependent flaring on dust growth in the 
future work. 

Figure[2]shows the total ionization rate ( as well as its three 
components (Ccr, Cxr, and Cra) as a function of r and z- X- 
ray ionization is dominant at outer radii (r > 4AU), while 
radionuclide ionization dominates at inner radii (r < lAU). 
Cosmic -ray ionization is important in outer (r > 2AU) and 
low-altitude (z < H) regions. 

3.1.3. Dust growth model 

Based on the results of recent laboratory and computer sim- 
ulations mentioned in § 1 , we model the dust growth in proto- 
planetary disks as follows. We start with monodisperse, non- 
aggregated dust grains (monomers) with radius ao. The mass 



niQ and number density ndo of monomers are then written as 
mo = (47r/3)poflo "do = Pd/'Wo- The dust is assumed to 
grow into an ensemble of quasi-monodisperse, fractal aggre- 
gates with typical monomer number and fractal dimension 
D ~ 2. Under this assumption, we may regard as the la- 
bel of dust evolutionary stages. Theoretically, this type of 
growth is the best modeled by th e so-called ba llistic cluster- 
cluster aggregation (BCCA; e.g.. iMeakinll 1 99 ID . For this rea- 
son, we shall refer to the quasi-monodisperse fractal growth 
as the "BCCA growth."' Note that the number density n^ of 
aggregates is inversely proportional to A^, since n^ = Nn^ is 
conserved. The fractal growth continues until colliding ag- 
gregates become energetic enough to c ompress each other. 
According to the microscopic model of iDominik & TielensI 
(1997), the critical kinetic energy for the onset of collisional 
compression is given by Lkin ~ ^roU, where 

2 /- 



^5.9x 10" 



f^)f^)erg (50) 
V2A/V0.1um/ ^ 



lOOerg/cm^/ V 2A ^ Vo.l^m 



is the energy needed to roll a monomer on another monomer 
in contact by 90°, 7 is the surface adhesion energy for the 
two monomers, and ^crit is the critical tangential displace- 
ment for starting the rollin g. For icy mono mers, 7 is es- 
timated as 7 « lOOg/cm^ dlsraelachvilil 11992) but a realis- 
tic value of ^crit is unknown. For a conservative estimation, 
we assume the minimum displacement £crit = 2A anticipated 
by the theory ( Dominik & Tielens 1997), which makes our 
aggregates the most easily compressed. The assumed value 
of iiroii is not much different from the experimental value 
for rocky (SiOa) monomers, Erow « 1.3 x 10~^(ao/0.1/im)erg 
dHeim et alj|1999 ). so the duration of the fractal growth stage 
is insensitive to our choice of dust material. We restrict our 
calculation to an early growth stage where the relative kinetic 
energy E]^in does not exceed the critical rolling energy E^on- 

The radius a of a fractal aggregate is approximately given 
by 

awaoA^'/^. (51) 

A classical, compac t aggregate ha s D = 3, while a BCCA 
cluster has D w 1.9 (IMeakinll 19911) . We adopt D = 2 in the 
following calculation. The projected cross section cr of an 
aggregate is simply set to a ~ na^ « aoN^^'^ w (7oN, where 
fJo = TTfly is the geometrical cross section of a monomer. This 
assum ption is consistent with the fact that a e x for D < 
2 (M eakin & DonnI Il988t iMeakin et alj [19891: iMinato et alj 
i2006) . Note that the quantity an^ is independent of A^, i.e., 
conserved for the BCCA growth. 

Assuming a quasi-monodisperse size distribution, the ki- 
netic energy ( [39] l is written as 



^kin ' 



\m(Au)\ 



(52) 



' Exactly speaking, the BCCA (i.e., collision between identical aggre- 
gates; see Meakin 1991) can only occur when the relative velocity of ag- 
gregates is induced by Brownian motion. When the relative velocity is in- 
duced by differential sedimentation, any identical aggregates cannot collide 
with each other, and therefore a BCCA cluster in its original sense cannot be 
created. On the other hand, a laboratory experiment by Blum et al. ( 1998|) 
shows that the outcome of sedimentation-driven coagulation is an ensemble 
of quasi-monodisperse, fractal aggregates with D K, \ as is for Brownian- 
motion-driven growth (e.g., Kempf et al. 1999). This may be explained by the 
fact that the dominant growth mode in differential sedimentation is the col- 
lision between similar aggregates iTanaka et al. 2005). For this reason, this 
study treats the growth by differential sedimentation as the "BCCA growth." 
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Fig. 2. — Total ionization rate ^ {thick gray curves) at different di.sk radii (a) and and altitudes (b). Here H denotes the scale height of the disk. The solid, 
dashed, and dotted black curves represent the contribution from Galactic cosmic rays, stellar X-rays, and radionuclides, respectively. 



where we have used m k, m/2. In a protoplanetary disk, rel- 
ative motion of aggregates is induced by Brownian motion, 
sedimentation toward the midplane of the disk, and turbu- 
lence. We therefore write the relative velocity Am as 



Am ~ V( AMBi-own)2 + (AMsed)^ + ( AMturb)^ , (53) 

where AMBrown, AMsed, and Awturb are the relative speed in- 
duced by the Brownian motion, differential sedimentation, 
and turbulence, respectively. 

The mean relative speed of the Brownian motion is given 
by 



AwBro 



16/tBr 



(54) 



■Km V 7rm 

In fact, the relative speed of the Brownian motion fluc- 
tuates according to the Maxwell distribution, and aggre- 
gates have a chance to get a relative kinetic energy E much 
larger than the thermal energy ^ k^T with a probability cx 
E^^^expi—E/k^T). However, as we see later, the effect of the 
thermal velocity fluctuation is insignificant, since the electro- 
static energy can go up to IO^^b^. 
The relative speed induced by the differential sedimentation 

is 



AMsed = f^K^Afs 



stop J 



(55) 

where fstop is the stopping time of an aggregate. For aggre- 



gates smaller than the mean free path ig of gas particles, 
is given by Epstein's law 

3 m 
4pgMg a 



'stop 



(56) 



where Mg = ^ 8A;b T / Tr/iwe is the mean thermal speed of gas 
particles. The mean free path in our disk model is calculated 
to be ig ^ l(r/l AU)"'''*cm, which is much larger than a typ- 
ical size of aggregate which we are interested in. For this 
reason, we always use Epstein's law ( |56] | in the following cal- 
culation. Also, we replace Afstop with its maximum value fstop, 
which leads to the most conservative evaluation of dust charg- 
ing effect. 

The relative velocity induced by turbulence is given by 



AMturb ~ Af, 

'small 



stop. 



(57) 



where Msmaii and fsmau are the characteristic velocity and 
turnover time of the smallest turbulent eddies , respectively 
jWeidenschilUnd [T984t lOrmel & Cuzzil l2007h . This ex- 
pression is valid for aggregates with stopping times fstop 



much smaller than f small- The velocity and turnover time 
of the largest eddies, Miarge and fiarge, are related to Msmaii 

and fsmaii by Miarge = Re'/'^Msmaii and fiarge ~ Re'^^fsmaii, wherc 
Re = i/tuib I Vmo\ is the Reynolds number. The molecular 
viscosity v^^x is written as fmoi = 0.5Mg/ngCrn,oi, where 
Omo\ = 2 X lO'^^cm^ is the molecular collision cross section 
jChapman & Cowiina,1970) . We express the turbulence vis- 
cosity I/turb = wLge^large aS f^turb = "turbcJl^K ' whcrC Qfturb is 

the a-parameter. The turnover time of the largest eddies is 
taken to be fiarge ~ ' ^'^'^ '■^^^ their velocity is given by 
Miarge ~ v^OtorbCj. We considcr aturb = 0, 10"^, 10"^, and 10"^ 
in this study. Again, we replace Afstop in equation (|57] | with 

^stop. 

If both of colliding aggregates have the mean charge (Z), 
the electrostatic energy (|40] | is written as 



2fl 



2A2' 



(58) 



where F is the master parameter defined in equation (|23] ). 
Equation ( fSST l overestimates a true repulsion energy when one 
(or both) of the aggregates has a positive charge Z > 0, or a 
negative charge -Z > smaller than the average value -(Z). 
To account for the dispersion of the charge distribution, we 
introduce the "three-sigma" electrostatic energy 



_ (Z)((Z) + 3(AZ^)'/V 
La 



(59) 



where the dispersion [l^Z^Y/^ is calculated from equation 
( l24l i. /iei.scr represents the electrostatic energy between two 
aggregates with charges (Z) and (Z) +3(AZ^)'''^. The prob- 
ability that a collision involves the electrostatic energy larger 
than -£'ei,3CT is as small as "three sigma" 10"^). Note that 
fieiSCT is always smaller than E^\ since (Z) is always nega- 
tive. Also, /iei.30- becomes negative if -(Z) is smaller than 
3(AZ^)'''^. This actually happens when the aggregate size is 
sufficiently small (see fig.[3]2 below). 

We compute the charge state at each evolutionary stage 
from equation ( l34l l with a w a{N), a w u{N). The sticking 
coefficients a re estimated by ph onon theory to be Si w 1 and 
ie = 0.1 . . . 1 (IUmebayashill983h . We adopt Si = 1 and = 0.3 
in this study. The results obtained in this section are insensi- 
tive to the choice of s^. as long as 0.1 < ie ^ 1- We do not use 
the original equation ( [30] l because the gas-phase ionization is 
negligible in the present case. The gas-phase recombination 
rate is typically /? 10"'^ ■ "^cm^/s. Using this value, we can 
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estimate g{r) as 



MiSeMe(0'«d)^ 
/c 



\io-y \5AuJ Vbok^ VlO-iVs^' 

(60) 

independently of A^. Therefore, equation (|34] | is valid unless 
dust is depleted and /jg decreases by many orders of magni- 
tude. 

To confirm that our semianalytical calculation does work 
well, we have also performed fully numerical calculations 
including multi-component ions. In t he numerical calcu- 
lations, a simple reaction model by lUmebayashi & Nakanol 
dTgSO.) is adopted. This reaction model involves five light 
ions (H"'",H2,H3,He''",C"'"), heavy molecular ions (m"*"), metal 
ions (M"""), free electrons, and charged dust aggregates. Heavy 
molecular ions and metal ions are represented by HCO""" and 
Mg"^, respectively. We adopt the same values of the rate co- 
efficients /3'^*\ P ^''''^ and the neutr al gas abundances n^''^ /n„ 



as those used bv lSano et al.l (l200d). As seen below, the dom- 
inant ions are metal ions, which is essentially due to the fast 
charge transfer from heavy molecules to metal atoms. With 
this fact, we set the average ion mass nii to be the mass of 
Mg"^ (mi = 24mu) in semianalytical calculations. The numeri- 
cal solutions are obtained from equations ([TJ-llSll. We remark 
that the full numerical calculation is far more time-consuming 
than the semianalytical one. 

3.2. Results 
3.2. 1 . The fiducial case: amrb = 0, oq = 0. 1 p.m 

Here we show the result for the case atmb = (i.e., laminar 
disk) and ao = 0.1 /im as a fiducial example. The results for 
different values of ao and aturb are examined in §3.2.2 and 
§3.2.3, respectively. 

Figure |3] shows the evolution of dust charge state and gas 
ionization state at r = 5 AU, z = H. Here each evolutionary 
stage is labeled by the number of constituent monomers in 
an aggregate, A^. The evolution of the mean (Z) and dis- 
persion (AZ-) of the dust charge distribution as well as the 
master parameter F is shown in figure [3]3. We find that F in- 
creases with and reaches to the maximum value Fmax = 2.8 
at N ^ 10^. This means that the gas-dust mixture is an ion- 
dust plasma (9 <C Fmax) at the initial stage, and evolves into 
an ion-electron plasma (Q 3> Fmax) as the dust grows. This is 
expected from the analysis in the last section: in the BCCA 
growth, 9 cx l/{aan^) is proportional to A^'/^, and thus in- 
creases with the growth. We have confirmed that the above 
value of Fmax is consistent with the solution to equation (l35T l. 

It is useful to introduce the critical monomer number A^ = 
A^max at which the transition from the ion-dust plasma regime 
to the ion-electron plasma regime occurs. As explained in 
§2.3, this value can be estimated by setting 9 « Fmax, or 



^^iFmax 



(61) 



Substituting w n^o/N, a w aoN, and a sa aoN^I^ into this 
equation, we obtain the critical size for the transition 



/ 2 T 

/ (T()floMinj,)i max 



1/2 



10 VI 



0.014 



5AU 



c 



130K/ VlO-iVs^ 



^0.1/im 
or equivalently, 



\0.\um) V1.4g/cmV 



(62) 



iaoN. 



1/2 
max 



^^^■^^(0.014) (sAu) (i30k)(iO-iVs) 
xf- 

vo. 



-1/2 



cm, 



(63) 



AfimJ V1.4g/cm3/ 
where we have explicitly expressed the dependence on /s. 

The mean (negative) charge -(Z) is proportional to A^ in 
the ion-dust regime (A^ <^ A^max), and is proportional to A^'/^ 
in the ion-electron regime (A^ ^ A'max)- This is easily under- 
stood if one recalls that -(Z) oc Fa oc FA^'/^. For A^ < A^max, 
F w 9 is proportional to A^'/^, so -(Z) cx A^. For A' > A'max, 
F approaches a constant, and thus -(Z) cx A^'/^. On the other 
hand, the dispersion (AZ^)'/^ is found to be nearly propor- 
tional to A^'/"^, which is because (AZ^)'/^ (a/A)'/^ and 
a cx A^'/^. It is important to notice here that the relative width 
I (AZ^)'''^/(Z) I of the charge distribution becomes sharper 
and sharper as the dust grows. 

The transition of the plasma state is better illustrated by fig- 
ure [3f'- This figure shows the number densities of ions and 
electrons, rii and «e, as well as the net dust charge density 
{Z)ni. For A^ ^ A^max, dust is the dominant carrier of nega- 
tive charges, as expected for the ion-dust plasma state. The 
free electron density is smaller than that of ions by a factor 
of (l/.Se)('«e/'Mi)'''^ 10"^ (see eq. |l38l)- As the aggregates 
grow and A^ reaches the critical number A^max, the dust charge 
density -(Z)nd begins to decrease and begins to increase. 
Finally, atN ^ 10^, free electrons become the dominant nega- 
tive charge carrier, and the ion-electron plasma state («i « Wg) 
is established. 

Interestingly, the abundances of charged species, n-, and n^, 
are nearly constant for both limits of ion-dust and ion-electron 
plasma regimes. This result is in contrast to that of previous 
studi es based on the classical, compact (D = 3) growth model 
(e.g., ISano et aT 2000; Wardle 2007) in which rii and /ig in- 
crease as the dust grows. This difference is attributed to the 
fact that the net projected area (T«d of dust aggregates is kept 
nearly constant for D <2, while it decreases for D = 3. Us- 
ing the constancy atii = croHdo, equations ( |32] | and ( |33] ) can be 
rewritten as rii = "io /(I + F) and «e = nio(ieMe/ Mi) expF, where 



MiCTOndO 

\0.0UJ V5AU/ VlO-'Vs/ 



\0.lnmJ Vl 



Po 



, / M 4 / ii (64) 
l/im/ \1.4g/cm-^/ 

is the abundance of ions for A^ ^ A^max- For both plasma lim- 
its, the factors 1 -l-F and expF are approximately constant, so 
both Mi and approach constant values. Physically speaking, 
the constancy of crwd for D <2 means that all the monomers 
in a fractal aggregate with D <2 are exposed to outer space, 
and are thus capable to capture free electrons and ions. 

Figure[3]also shows the result of full numerical calculations 
using the simplified ion-reaction scheme. We find an excel- 
lent agreement between the semianalytical and full numerical 
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Fig. 3. — (a) The evolution of the dust charge distribution for the "BCCA" (i.e., quasi-monodisperse, D 2) dust growth. The solid, dashed, and dotted curves 
represent the average (Z), dispersion (AZ^)'/^, and master parameter F of the charge distribution calculated using the semianalytical method (eqs. l23l . 1241 . and 
1341 ). The filled circles and squares denote (Z) and (AZ^)'^~ obtained from full-numerical calculations using the Umebayashi & Nakano model (see §3.1.3). (h) 
The evolution of the gas ionization state and the dust charge state for the BCCA growth, ni and are the number densities of ions and electrons in the gas phase, 
and (Z)fid is the net dust charge density. nf^+ and n„,+ are the number densities of metal ions and molecular ions. 
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Fig. 4. — The relative kinetic energy fiy,, (eq. 1521 ; solid black curve) 
versus the electrostatic energy (eq. [(58)]; solid gray curve) for the 
BCCA growth as a function of the monomer number N. The disk position 
is (r,z) — (5AU,i/), and the turbulence parameter and the monomer size are 
set to Qturb = (i.e., laminar) and cio = O.l/^m. The collision cross secti on fo r 
an aggregate pair with and vanishes when E^x > ZJ^in (see eq. 1421 ). 
meaning that the BCCA growth "freezes out" at the size indicated by the cross 
(x) symbol. Th e das hed gray curve shows the "thi'ee-sigma" electrostatic en- 
ergy Ziel.So- (sq- 1591 ), representing the effect of charge fluctuation. The star 
(*) symbol indica tes the size at which fiy^ reaches the critical rolling-friction 
energy iJrou (eq. 1501 ; dotted black curve). Above this critical energy, coUi- 
sional compression of aggregates becomes effective. 

calculations. It is clear that our semianalytical approach is not 
only efficient but also accurate. The most abundant ions are 
metal ions for all stages of dust evolution. Molecular ions, the 
second most abundant ones, are an order of magnitude fewer 
than metal ions. 

Now we examine the growth condition. Figure |4] shows 
the kinetic energy /Skin and electrostatic energy /Sei for collid- 
ing aggregates at each evolutionary stages. For < 10^ (a < 
30fim), the thermal (Brownian) motion dominates the rela- 



tive velocity of the aggregates, so Ey^n is kept constant « k^T. 
Fore > 10'', the vertical sedimentation dominates the aggre- 
gate motion and fikin increases with A^. On the other hand, ^ei 
always grows with A^, in proportional to N^^- for A^ < A^max 
and to A^'/^ for A^ > A^max- This is explained from the fact 
that I (Z) I oc A^ for the ion-dust plasma regime (A^ < A^max) 
and |(Z)| oc fl cx 

A^i/2 for the ion-electron plasma regime 
(A^ > A^max)- As a result, the growth condition (02]) breaks 
down at N ^ 10''-^(fl w 6/im). This means that the collision 
between aggregates with average charge (Z) becomes impos- 
sible at this stage. Note that the repulsion energy Ef.\ reaches 
10 times the thermal energy ^ k^T as early as N ^ 10"*, and 
finally goes up to IO^^b?^ at the onset of collisional compres- 
sion. It is evident that the thermal fluctuation of the kinetic 
energy cannot help the aggregates to grow beyond A^ 10"*. 

One may expect that the fluctuation of aggregate charge 
could help the growth. To see this effect, we overplot in 
figure m the "three-sigma" electrostatic energy iSei.Scr defined 
in equation ( |59l ). We see that /iei.Scr quickly converges to 
/iei and finally exceeds iskin at N 10*. This is expected 
from figure [3}2: the relative width \{Z) / {AZ^)^/^\ of the the 
charge distribution becomes narrower and narrower as the 
dust grows. Therefore, the result that the aggregates cannot 
grow beyond A^ ^ 10"* is preserved even if the charge fluctua- 
tion is taken into account.^ All the above facts suggest that the 
BCCA growth of dust aggregates at this disk position "freezes 
out" at size A^ ^ I0'^(a ~ lO^m). 

It is interesting to compare this critical size for the freeze- 
out with that for the first compression, i.e., the size at which 
the growth mode changes from the BCCA to the growth in- 
volving collisional compression. We overplot the critical 

It is also found that E^\ j^ is negative for smaller sizes, N < 10"". This 
means that colliding aggregates can possess opposite charges with a proba- 
bility larger than "three-sigma". This is because the average negative charge 
-(Z) in this stage is smaller than 3(AZ-)'/^, as seen in figure[3j(. However, 
the attraction energy -Scl.So- is insignificant: at most half of the kinetic energy 
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rolling energy iiioii (eq. ISOl ) in figure|4] Comparing this crit- 
ical energy with Ekin, we find that the collisional compression 
begins at size « 10"'(fl w 1cm), i.e., many orders of mag- 
nitude smaller than the above critical freeze-out size. This 
illustrates how fast the charging of aggregates begins to affect 
their collisional growth. 

Figure |5] compares the energy ratio iSei/^kin for different 
disk positions. We find that the barrier against the BCCA 
growth appears irrespectively of the disk radius r. By con- 
trast, the growth barrier vanishes at z > 3// because the sed- 
imentation velocity at the high altitudes is large enough for 
aggregates to overcome the barrier 

To summarize, the freeze-out of the BCCA growth is very 
likely to occur in this fiducial model, except at high altitudes 
over the midplane. 

As seen in figure |4] the energy ratio /iei/^^kin takes its max- 
imum at « A^max- With this fact, we can roughly estimate 
the maximum value (-£'ei/£kin)max of the energy ratio at z w // 
as follows. We assume that the relative motion of aggregates 
by w Ntrna is dominated by vertical sedimentation, as is for 
r = 5AU. Then, substituting Am w AMsed into equation ( |52] | 
and using equations ( ISSl l and ( |62] |. we obtain 




at z w //. Notably, (^ei/^^kiiOmax is explicitly independent of 
both r and T . This means that the "height" of the growth 
barrier is insensitive to the temperature profile. Equation ( l65T l 
does not hold beyond r « 5AU, since the maximum energy 
ratio appears in the Brownian motion regime. This equation 
is nevertheless useful because it allows us a rough estimation 
on how the growth barrier depend on the model parameters. 
For example, for fixed C, qq, and po, the the growth barrier 
is more serious if the disk is more massive (/s > 1) or more 
depleted of dust (/dg < 10"^). 

3.2.2. Effect of monomer size 

The actual size of dust monomers in protoplanetary disks is 
unknown. Infrared observations of interstellar medium sug- 
gest the size distribution of interstellar grains ranges from 
» 0.0 05Mm to w 0.25^m (MRN disti-ibution; iMathis et all 
Il977h . If interstellar grains are not aggregates but monomers, 
the typical monomer size in protoplanetary disks will fall 
within the range 0.001/im < a < 1/im. Figure |6] compares 
the energy ratio /iei/^^kin for different monomer sizes = 
0.01, 0.1, and 1/xm. We find that A^^ax oc a^^ (flmax oc a^^) 
and (iiei/^kiiOmax OC Qq'' ^s cxpccted from equations (|62] | and 
i65[ . This delay (measured in the "degree of growth" A^) is 
attributed to the larger cross section of BCCA clusters com- 
posed of smaller monomers. Such a larger cross section 
causes a quick depletion of free electrons in the gas phase, 
resulting in the delay of transition from the ion-dust plasma 
regime to the ion-electron regime. In addition, the larger 
cross section produces a stronger coupling to the gas, and 
in turn a slower increase of the kinetic energy. We find that 
the growth condition (|42] | breaks down much before the on- 
set of collisional compression for all aq < Ifim. Therefore, 
the "freeze-out" of the BCCA growth is not prevented even 
if the monomer size in protoplanetary disks is assumed to the 
maximum value inferred by the MRN distribution. 



3.2.3. Effect of turbulence 

Here we examine how strong turbulence is needed to re- 
move the growth barrier Figure [T] shows the evolution of en- 
ergy ratio Eei/^kin for three turbulent cases aturb = 10"^, lO"-', 
and 10"^. We see that turbulence of aturb ^ 10""^ does not 
affect the energy ratio for any size A^. This is not surpris- 
ing because both AMtm-b and Awsed scale with Afstop- The 
relative velocity Awturb induced by turbulence is estimated 

as AMturb ~ AfstopRe'/'^Miarge/flarge ~ At ^top^S^ ^'^ aHlnl^H . At 

H, the Reynolds number Re is of order ~ amrbSgO-moi / 
so AMtui-b is written as 

A«,urb - (^^) '^'a,Y>f,opf]ii/ (66) 



On the other hand, the relative velocity by differential sedi- 
mentation is AMsed ~ At^iopft^H, so we find 




Therefore, the effect of turbulence on the aggregate collision 
is negligible for all A^ as long as a^h ^ 10 . 

For Qfturb ^ 10"^, we find that the growth barrier is entirely 
removed. This suggests that relatively strong (atuit ^ 10"^) 
turbulence is a key ingredient for early stages of dust coagu- 
lation. In §4.3, we discuss this topic in more detail. 

4. DISCUSSION 
4.1. Validity of the charge equilibrium 

In this study, we have assumed that the charge reactions are 
much faster than the dust-dust collisions. Now we show that 
this assumption is actually valid for evolutionary stages which 
we are interested in. 

The typical time scale for the system to relax to an 
ionization-recombination equilibrium can be measured by the 
average time fcoii.d needed for an aggregate to collide with an 
ion, 

-1 C'^e 

fcoU i ~ «iC^di«i « , (68) 

«d 

where we have used that (^n„ k, Uiannin^ since the gas- 
phase recombination is negligible in the presence of dust (see 
§3.1.3). On the other hand, the mean collision time fcoU.d be- 
tween aggregates is written as 

Qn,d ~ CTddHdAM ^ (Tn^Au. (69) 
Therefore, the ratio of these time scales is estimated as 

fcoU.d C»g .-_ f-2f"^ V '"b ^ (70) 

fcoU.i an^Au '^°\mgJ Am crSg J^k 

If the dust velocity is dominated by the Brownian motion, as 
is for small aggregates. Am ^ c^^jmjm and 




where we have used aq ^ O.l/im and po ~ Ig/cm-*. We 
find that dust coagulation can be safely neglected if r > 2AU 
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Fig. 5. — The ratios E^xjEy^^ between the electrostatic and kinetic energies (solid curves) at various disk radii r (a) and altitudes z (h) as functions of A?. The 
o-parameter and the monomer radius are set to Qtmb = " ^nd oq = 0.1/im. The cross (x) symbols indicate the size at which the growth condition i42i breaks 
down for aggregates with mean charge (Z). The dashed curves show the ratio fid Str/^^kin for the "three-sigma" energy, which represents the effect of charge 
dispersion. The star (*) symbols indicate the sizes at which fi^j, reaches the critical rolling-friction energy Smll- 
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Fig. 6. — The energy ratios i?cl/£^kin (solid curves) at r = 5AU and z = H 
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out" sizes at which the growth condition i42) breaks down. The dashed 
curves show fielScr/Skin- The star (*) symbols indicate the sizes at which 
reaches the critical rolhng-friction energy E^^w- 



or N ^ 1. This is true even if the motion of dust aggre- 
gates is dominated by vertical sedimentation, since Awsed ~ 
Cs{m/Y,ga), and thus 




?coll,d 



fcolli 



10 V 



m C 

/dg 



■''^^ mg riK 



10-2 



lO-'Vs 



(2Au) 



3/2 



(72) 



Thus, it is concluded that the growth of dust aggregates can 
be neglected if r > 2 AU ovN:s> I. 



Fig. 7. — The energy ratios £cl/£kin (solid curves) at r = 5AU and z = 
H for different turbulence strengths citurb- The cros s (x ) symbols indicate 
the "freeze-out" sizes at which the growth condition )42t breaks down. The 
dashed curves show £el.3o-/£^kin- The star (*) symbols indicate the sizes at 
which iJitii, reaches the critical rolling-friction energy E^^w- 



4.2. Internal electrostatic force 

In the last section, we have implicitly assumed that charged 
aggregates can stick to each other as long as the collision con- 
dition ( |42] | is satisfied. One might wonder if the collided ag- 
gregates are pulled off from each other by the electrostatic 
repulsion. In fact, this repulsion is much weaker than the at- 
traction (due to van der Waals force) between two monomers 
in contact. The electrostatic repulsion force Fei.int acting be- 
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tween two collided aggregates is estimated as 



{Zef 



TV 



^elint ' 



< 10" 



dyn, 



(73) 



fl2 A2 V130K. 
where we have used that F < r^ax ^ 3. Note that max- 
imum value of int is independent of the aggregate size 
a. On the other hand, the critical force needed to sep- 
arate two monomers m contact is (iJohnson et alJ 119711: 
iDominik & Tielenslll997h 



10" 



2 " VlOVanV VO.l^my 

where 7 is the surface adhesive energy mentioned in the last 
section. Thus, the electrostatic force inside an aggregate is 
negligibly weak compared to the contact force between two 
constituent monomers. 

4.3. Dust growth in strong turbulence 

As seen in §3.2.3, the electrostatic barrier against the fractal 
growth will be removed if considerably strong (aturb }t 10"^) 
turbulence is present. We here discuss whether such turbu- 
lence is likely to occur in protoplanetary disks, and what 
would happen after the dust overcome the electrostatic growth 
barrier. 

The most robust mechanism for disk turbulence is MRl ^. 
MRl-driven turb ulence will achie ve amib ^ 10"^ in its satu- 
rated state (e.g. ISano et aT]|1998l) . Therefore, fractal aggre- 
gates will be able to overcome the electric barrier in MRI- 
active regions. ISano et all (1200 0) calculated the active region 
using the MMSN model and found that MRI will be active 
only at outer (r > 20AU) disk radii or high (|z| > 2H) altitudes 
if the dust size is 0.1 /im. The size of the active region does 
not vary even if the dust grows since, as seen in §3.2.1, the 
ionization fraction is kept nearly constant as long as the dust 
undergoes the fractal growth. Therefore, the region in which 
the fractal aggregates can overcome the electrostatic barrier is 
limited to outer disk radii and high altitudes. 

A more serious problem is that such strong turbulence 
causes another kind of growth barrier, i.e., catastrophic frag- 
mentation of colliding aggregates. In turbulent regions, large 
aggregates with fstop ^ 1 /^^k have the maximum collisional 
velocity of order Miaige ~ ^/atm-bCj, which amounts to more 
than lOOm/s for atuib ^ 10"^ . On the other ha nd, as shown 
by recent A^-body simulations (IWada et al.l2008h . catastrophic 
fragmentation will take place for relative velocity Am > 
30m/s. Therefore, it is very likely that strong turbulence de- 
stroys the aggregates and consequently prevents further dust 
growth. This idea is supported by a recent statistical study 
iSrauer et al.1l2008h . 

Thus, the combination of electric repulsion and collisional 
fragmentation might strictly limit the dust growth and sub- 
sequent planetesimal formation in protoplanetary disks. It 
is important to think of a possibility that dust evolution will 
continue in some way even if the turbulence is weak and the 
quasi-monodisperse fractal growth does freeze out. This is 
the topic of the next subsection. 

4.4. A possible scenario to overcome the electric growth 

barrier 



' It is unknown whether any mechanism other than MRI can drive and 
sustain turbulence with cttuib ^ 10^ in ths early stage of dust evolution. For 
exa mple, convective instability may operate in this stage ( Lin & Papaloizou 
M982^, but it is unlikely to sustain such strong turbulence iStone & Balbusi 



In §3, we have ignored the size distribution of aggre- 
gates. In fact, there may exist some aggregates considerably 
larger than average-sized ones. In the following, we consider 
whether such large aggregates can continue to grow even if 
the growth of average-sized ones has frozen out. 

Let us consider a small population of irregularly large ag- 
gregates (referred to as "test aggregates") growing with a large 
population of standard {D ^ 2) fractal aggregates ("field ag- 
gregates"). Under this assumption, the kinetic energy of rela- 
tive motion between test and field aggregates is written as 



1 niimf 



2 nil + ni{ 



(AMtf)2 « -mf(AMtf)2 =-' " ^ 



2 V Awff / 



£'kin,ff, 

(75) 

where the subscripts 't' and 'f respectively represent the test 
and field aggregates, and we have used the assumption ^ 
mf. iikin.if and Aug are the kinetic energy of relative motion 
and the relative velocity between two field aggregates, and are 
thus equivalent to fikin and Am in §3. On the other hand, the 
electrostatic energy between test and field aggregates is 



^el.tf = 



aiafT-e 



2„2 



r^flf 1 



(76) 



where ZSeiff is equivalent to E^i in §3, and we have used that 
flt S> flf. Thus, the energy ratio -Eei.tf/^^kin.tf is written as 

£'el,ft £'el.ff /AMff\2 



^kin.ft ^kinjf ^ AMtf / 



(77) 



Now we assume that the growth of field aggregates has 
frozen out due to the electric barrier, i.e., /ikin.ff = £ei,ff- At this 
stage, the condition for the collision between test and field 
aggregates, £'kin,tf > £ei,tf, reduces to a simple inequality 



AMtf > AMff. 



(78) 



We readily notice that Brownian motion does not sat- 
isfy this condition since AMtf ~ y/ 8^b T / Trnif and AMff w 

^l6ksT /Trnif. The remained possibilities are the differen- 
tial sedimentation and turbulent-driven motion. In both cases, 
the relative velocity Am is proportional to Afstop, or A(m/a). 
It is very important to notice here that the condition ( fTsT l is 
not safely satisfied as long as the test aggregate is as fluffy 
as field aggregates, i.e., mjat is comparable to nif/af. Hence, 
the condition ( iTSl l will be safely satisfied only if the test aggre- 
gate is more compact and has larger m/a than the field aggre- 
gates. Moreover, collision with a smaller aggregate generally 
tends to increase nit/ at, allowing the test aggregate the next 
collision. Therefore, if there exists an aggregate that is large 
and compact, it will be able to continue growing by sweeping 
up smaller "frozen" aggregates. 

The above consideration suggests that the freeze-out of the 
quasi-monodisperse fractal growth may not mean the termi- 
nation of dust evolution. Rather, it may be the beginning of 
bimodal growth in which only a small fraction of aggregates 
can grow larger and larger while the rest remain frozen. We 
plan to examine this possibility in more detail in the future 
studies. In any case, we expect that the effect of dust charg- 
ing should qualitatively modify the current scenario of dust 
growth in protoplanetary disks. 

5. SUMMARY 

In this study, we have investigated the electric charging of 
dust aggregates and its effect on collisional dust growth in 
protoplanetary disks. We have found that the conditions for 
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ionization-recombination equilibrium are reduced to a single 
equation (eq. |[30l ). Just by solving this equation numerically, 
the dust charge state and gas ionization state can be analyti- 
cally computed for an arbitrary ensemble of aggregates in a 
self-consistent way. It is also confirmed that our semianalyti- 
cal method reproduces the results of a previously used, more 
complicated numerical method (§3.2.1, fig. 3). This formal- 
ism thus provides a fast charge-state solver that will allow a 
coupled simulation of MRI-driven turbulence and dust coag- 
ulation. 

As an application, we have explored the effect of electro- 
static charging on an early stage of dust coagulation in proto- 
planetary disks. We considered the quasi-monodisperse frac- 
tal growth with the fractal dimension D ~ 2 as suggested 
by previous laboratory experiments and A^-body simulations 
(jBlum 2004; Dominik et al. 2007). Our findings are summa- 
rized as follows: 

1 . For a wide range of model parameters, the effective cross 
section for the mutual collision of aggregates is quickly sup- 
pressed as the fractal growth proceeds and finally vanishes 
at a certain aggregate size (§§3.2.1, 3.2.2). This is due to 
the strong electrostatic repulsion between aggregates charg- 
ing negatively on average, and happens much before the col- 
lisional compression of aggregates becomes effective. Both 
the charge fluctuation and the thermal velocity fluctuation do 
not help the aggregates to overcome the growth barrier With- 
out strong turbulence, the quasi-monodisperse fractal growth 
is very likely to "freeze out" on its way to the subsequent 
growth stage. 

2. Strong (atuib ^ 10"^) turbulence will help the aggregates 
to overcome the above growth barrier (§3.2.3). However, such 
turbulence is likely to occur only in MRI-active regions, i.e.. 



at outer disk radii or high altitudes (§4.3). Furthermore, it 
will cause another serious problem — the catastrophic disrup- 
tion of collided aggregates — in later stages. These facts sug- 
gest that the combination of electric repulsion and collisional 
disruption may strictly limit the collisional growth of dust ag- 
gregates in protoplanetary disks. 

3. The freeze-out of the fractal growth might be followed 
by bimodal growth in which only a small fraction of large ag- 
gregates can continue growing while a large fraction of small 
fractal aggregates remains frozen (§4.4). This could qualita- 
tively change the curre nt scenario of planete simal formation 
in protoplanetary disks (iDominik et al.ll2007h . We will exam- 
ine this possibility in more detail in forthcoming papers. 

Finally, we point out that the fractal (D < 2) dust growth 
tends to keep the ionization degree of the disk small due to 
the open nature of aggregates (§3.2.1, fig. 3b). This means 
that the magnetorotationally unstable region hardly expands 
until the collisional compression of the aggregates begins to 
work. This conclusion is in contrast t o that of previous studies 
(e.g., ' Sano et al.ll2000l: IWardlell2007 ) which claimed that the 
ionization degree increase as the aggregates grow. However, 
they assumed compact dust growth, which clearly contradicts 
recent laboratory experiments and A^-body simulations. Thus, 
the magnetorotational stability of protoplanetary disks must 
be reexamined taking into account that the fractal nature of 
dust aggregates. 
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APPENDIX 

DERIVATION OF THE CHARGE DISTRIBUTION 

Since the velocity of free electrons is much greater than that of ions, dust aggregates charge up negatively on average. Thus, 
let us assume Z < and use the expressions of (Tdi and crje valid for Z < 0. Then, equation (fTST l is written as 

/ Z\ fz+i\ 

SiUitli 1 nd(/,Z) = ieMeWeexp nd(I,Z+ 1), (Al) 



where t = a/X. 
We make the following assumption: 



' - <^^'>"^.. ,A2, 



(AZ2)l/2 T 

where (AZ^) is the variance of the charge state distribution, and e <C 1. We also assume that nd(/,Z) varies with a typical scale 
~ (AZ^)'/^. Under these assumptions, nd(^,Z+ 1) can be written as 

n,iI,Z+ 1) = «d(/,Z)+ ^(/,Z) + 0(e2). (A3) 
oZ 

Also, exp[(Z+ 1)/t] is written as e\p(Z / t) + O(e^) since r"' = (AZ^)"'/^ • (AZ^)'''^/t ^ Substituting them into equation 
dAU . we obtain a first-order differential equation for ni(I,Z), 

driA 
dZ 



^(I,Z) + W(I,Z)ni(I,Z)^0, (A4) 



where 

SgUgrie exp(Z/T) 



W(I,Z)=l--!-^ -f---. (A5) 



Equation ( IA4l i is accurate to terms of the first order in e. 

Let us denote the solution of W(/,Z) = by Zq and write Z = Zo + 6Z. Also, we make an additional approximation that \5Z\ ^ t. 
Expanding W(I,Z) in powers of \SZ\/t and using W(/,Zo) = 0, we have 

W(/,Z)«^^, (A6) 

l-Zo/T T 
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which is accurate to the first order in |i5Z|/r. Hence, equation ( IA4I ) is approximated by 

^(/,Z) + ^^«d(/,Z)«0, (A7) 

where 

(AZ^) ^ i^r. (A8) 

It is an easy task to show that the solution to ( |A7t is a Gaussian distribution (l22l i with average (Z) = Zq and variance (AZ^). 
Rewriting Zp and r using Z p = -F t and r = a/ A, we obtain equations ( |23] l and ( l24l i. We note that equation dASI l has been also 
obtained bv lDraine & SutinI (Il987h . but they did not show the derivation of this equation in their paper 

Now we justify the assumption (IA2l i. The first part of the assumption (AZ-)"'/^ ^ i^Z^) ^^^/t, which is equivalent to (AZ^) ~ 
r, is always satisfied since t/2 < (AZ^) < r (see eq. HASI ). The second part (AZ^)'''^/t ^ e <C 1 is also satisfied in typical 
protoplanetary disks if dust aggregates are not as small as constituent monomers Q.l^m), because 



(A9) 



It is noted that the approximation |(5Z| <C t used to obtain equation (IA6b is rewritten as \8Z\ <^ t^I'^{1\ZYI'^. Hence, this 
approximation is good as long as the region \5Z\ < (AZ)i/2 is considered. 
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